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ABSTRACT 


An  algebraic  reconstrnction  technique  (ART)  ia  described  for  the 
seisnic  tonography  of  velocities  fron  the  travel  tines  of  a  nnltiple 
offset  vertical  seisnic  profile.  ART  concentrates  on  the  production  of 
a  reconstructed  field  whose  projected  data  (travel  tines)  agree  with  the 
observed  data.  This  reconstructed  field  is  nodified  by  altering  the  data 
for  each  ray  such  that  when  this  data  is  back-projected,  the  new  inage 
agrees  with  the  original  data.  Because  the  paths  of  the  rays  nust  be 
known  to  calculate  the  expected  travel  tines,  the  problen  is  linearised 
by  using  raypath  approxinations  as  detemined  fron  either  a  constant  or 
a  linear  c(z)  velocity  nediun. 


Inaging  of  synthetic  data  revealed  that  the  orientation  of  the  ano- 
naly  affects  both  the  rate  of  convergence  and  the  resolution  of  the 
reconstructed  field.  Sone  snoothing  of  the  velocity  anonalies  occurred 
along  the  direction  of  the  rays. 


Noisy  data  sets  developed  problens  in  the  reconstructed  velocity 
field.  Huge  single  point  anonalies  appeared  along  the  nodel’s  edge  in 
the  reconstructed  inage./  Elinination  of  these  errant  anonalies  was 
necessary  to  obtain  a  reasonable  velocity  reconstruction.  Because  iso¬ 


& 


tropy  was  assuned.  the  algorithn  poorly  reconstructed  data  which  was 
collected  over  the  strongly  anisotropic  Pierre  Shale. 
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NOTATION 


Angle  forned  by  the  length  of  the  raypath. 


Velocity  at  the  tnrface. 

Velocity. 

Son  of  chord  lengths. 

Residual  error  time  for  chord  length. 

Residual  error  time  for  entire  ray. 

Angle  of  incidence. 

Subscript  -  from  the  i'th  source  to  the  j'th  receiver. 
Chord  length. 

Total  length  of  the  ray. 

Superscript  for  parameter  circle  around  point  a.n. 
Slope  of  the  linear  c(z)  velocity. 

Ray  parameter. 

Radius  of  any  parameter  circle. 

Radius  of  curvature  for  the  ray. 

Slowness  matrix. 

Slowness  correction  matrix. 

Residence  time. 

Calculated  travel  time  for  the  ray  i j . 

Observed  travel  time  for  the  ray  i j . 

Weighting  factor. 


v*j 


Offset  to  source,  depth  to  receiver. 
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1.  INTRODUCTION 

One  of  the  nany  goals  of  vertical  seisaic  profiling  (VSP)  is  to 
determine  the  lateral  variations  which  occur  away  from  the  borehole. 
These  lateral  changes  aay  reflect  porosity,  permeability,  or  facies 
changes  within  the  rock  unit.  Such  occurrences  greatly  influence  the 
migration  and  accumulation  of  petroleum.  Therefore,  methods  are  needed 
to  delineate  these  discrepancies  between  the  borehole  wall  and  the  aur- 
rounding  medium. 

This  investigation  was  geared  at  determining  velocity  inhomo¬ 

geneities  using  only  the  first  arrival  times  from  a  multiple  offset  VSP. 
The  travel  time  through  any  medium  is  the  integral  of  the  slowness 

(reciprocal  of  the  velocity)  along  the  raypath.  Unfortunately,  this  ray- 
path  itself  is  dependent  upon  the  velocities  within  the  medium.  Thus, 
as  posed,  the  problem  is  non-linear. 

Linearizing  this  problem  requires  that  the  raypaths  are  known, 

which,  in  turn,  necessitates  a  prior  knowledge  of  the  velocity  struc¬ 
ture.  However,  if  we  assume  some  general  velocity  trend  for  the  medium, 
we  can  approximate  the  paths  which  the  rays  take.  Then,  the  problem 

becomes  one  of  finding  the  velocity  perturbation  from  the  assumed  back¬ 
ground  velocity.  Furthermore,  we  assume  that  all  raypaths  are  not 


influenced  by  this  perturbation  velocity;  this  assumption  holds  true  for 
small  perturbations.  Since  the  rays  are  independent  of  the  velocity 
perturbation,  the  problem  becomes  linear. 

This  research  utilized  two  background  models: 

1)  Constant  velocity  medium  in  which  the  rays  form  straight  lines  from 
the  source  to  the  receiver. 

2)  Linear  c(z)  velocity  medium  where  the  raypaths  follow  arcs  along  a 
circle. 

The  residual  time  is  defined  as  the  observed  time  minus  the 
expected  time  as  calculated  from  the  velocity  model.  Fawcett  (1983) 
utilized  a  Radon  transform  method  for  his  seismic  tomography  of  the 
slowness  field  from  reflected  residual  travel  times.  Neumann  (1981) 
used  a  least-squares  inversion  of  residual  travel  times  in  his  study  on 
reflection  seismics.  Christoff erson  and  Bnsebye  (1979)  applied  a  least- 
squares  approach  on  P-wave  residual  times  in  a  three-dimensional  case. 
In  contrast.  McMechan  (1983)  and  Nason  (1981)  opted  for  an  algebraic 
reconstruction  technique  (ART) ,  I  have  opted  for  the  last  approach,  as 
well . 

In  comparision  to  a  least-squares  inversion,  ART  possesses  several 


advantages : 


1)  ART  programs  are  computationally  faater  and  they  are  eaaier  to  pro¬ 
gram. 

2)  Constraints  are  easily  incorporated  into  the  program  to  accommodate 
any  prior  knowledge  of  the  medium. 

3)  ART  can  be  easily  applied  to  any  source/receiver  geometry  without 
difficulties. 

ART  has  its  origins  within  the  medical  field.  Gordon.  Bender,  and 
Herman  (1970)  first  introduced  ART  for  the  reconstruction  of  images  from 
x-ray  pictures.  Later,  Herman,  Lent,  and  Rowland  (1973)  improved  the 
early  ART  algorithms  and  Gordon  (1974)  summarized  the  various  ART 
developments  within  the  medical  profession. 

Several  authors  worked  on  reconstruction  techniques  which  are 
applicable  to  the  more  general  reconstruction  problem  rather  than  the 
specific  x-ray  case.  Hersereau  and  Oppenheim  (1974)  and  Mueller,  lav eh, 
and  Wade  (1979)  applied  ART  to  density  reconstruction  problems  using  the 
Fourier  domain.  Horn  (1978)  developed  a  density  reconstruction  algo¬ 
rithm  for  any  arbitrary  ray  sampling  scheme. 

Geophysical  tomography  is  relatively  new.  Dines  and  Lytle  (1979) 
reconstructed  pictures  of  electromagnetic  properties  in  the  region 
between  a  pair  of  boreholes.  By  assuming  only  straight  raypaths, 
seismic  velocity  structures  were  obtained  by  both  Mason  (1981)  and 
McMechan  (1983).  Mason  concerned  himself  with  inverting  seismic  data 
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shot  between  two  boreholes,  while  McMechan  imaged  data  shot  from  a 
borehole  to  a  second  borehole  as  well  as  to  the  snrface. 

The  VSP  model  nsed  in  this  paper  consisted  of  ten  point  sources, 
which  were  equally  spaced  at  100  foot  intervals  along  the  snrface. 
Similarly,  a  100  foot  spacing  was  nsed  for  the  ten  receivers  placed  down 
the  borehole.  No  source  or  receiver  resided  at  the  top  of  the  borehole 
(see  Fig.  1).  For  a  coordinate  system,  whose  origin  rested  at  the  snr¬ 
face  expression  of  the  borehole,  the  x  variable  coincided  with  the  snr¬ 
face  offset  distance  and  the  z  variable  reflected  the  depth  below  the 
snrface. 

2.  ALGEBRAIC  RECONSTRUCTION  TECHNIQUE  -  ART 

The  ART  algorithm  nsed  here  was  adapted  from  Mason  (1981).  It 
involves  a  nine  step  process: 

1)  The  medium  is  modeled  by  casting  a  set  of  overlapping  circles  on  a 
uniform  grid. 

2)  The  total  length  of  ray  segments  abont  a  parameter  is  the  sum  of 
all  the  chord  lengths  formed  by  rays  which  pass  through  any  given 
parameter  circle. 

3)  The  travel  times  associated  with  each  of  these  chord  lengths  are 
determined  by  the  velocity  background  model  and  summed. 


4)  The  average  slowness  at  the  center  of  a  parameter  circle  is  the 
difference  of  2)  and  3).  These  outputs  form  a  matrix  indexed  by 


5)  Using  the  generated  slowness  matrix  from  4),  the  expected  travel 
times  are  calculated. 


6)  The  residual  time  is  the  observed  minus  the  calculated  travel  time. 


7)  This  residual  error  time  is  equally  distributed  along  the  entire 
length  of  the  ray. 


8)  A  slowness  correction  matrix  is  computed  in  the  manner  of  steps  2), 
3) ,  and  4) . 


9)  Steps  3)  through  8)  are  iterated. 


2.1  Constant 


As  stated  earlier,  this  investigation  used  two  different  velocity 
backgrounds  for  the  medium.  First  the  constant  velocity  case  will  be 
discussed  before  going  into  the  more  complicated  linear  c(x)  medium. 


Using  a  constant  velocity  reference,  the  raypaths  become  straight 
lines  from  the  i ' th  source  to  the  j'th  receiver.  Denote  the  length  of 
this  ray  as  L^.  Because  the  ray  follows  a  straight  course  from  (xj.O) 
to  (0,z.),  the  raypath  is  represented  by  the  equation: 


T 


Over  the  two-dimensional  medium,  ■  set  of  overlapping  circles  is 
drawn,  whose  centers  form  an  uniform  square  grid  (Fig.  2).  These  cir¬ 
cles  must  completely  cover  the  medium  such  that  every  point  within  the 
medium  is  contained  in  at  least  one  of  these  parameter  circles.  The 
average  slowness  within  any  circle  is  assumed  to  be  the  slowness  at  the 
origin  of  the  circle.  These  circles  serve  as  a  truncated  least-distance 
averaging  method  which  lies  at  the  basis  of  Horn's  algorithm  (1978). 
Also,  circles  prove  to  be  more  convenient  than  a  grid  of  squares. 

Sources  are  positioned  at  the  center  of  some,  or  all.  of  the  cir¬ 
cles  along  the  surface  and  receivers  are  located  at  the  origin  of  some, 
or  all,  of  the  circles  at  the  borehole  edge  of  the  model.  Clearly,  the 
source  or  receiver  spacing  defines  the  largest  possible  grid  size,  or 
equivalently,  the  density  of  the  parameter  circles.  Resolution  is  lost 
with  too  small  of  a  density  of  parameter  circles.  Conversely,  a  greater 
grid  density  will  increase  the  resolution  only  to  the  limit  of  the  data. 

Let  the  known  value  r  denote  the  radius  of  any  parameter  circle. 
Then  the  circle,  around  the  point  (x  ,z  ),  is  mathematically  expressed 

D  O 

as : 

(x-x  )*  +  (z-z  )*  =  r*.  (2) 

m  n 


If  the  ray  in  question  intersects  the  parameter  circle,  then  equa¬ 
tion  (1)  and  equation  (2)  vill  have  two  points  in  common.  A  quadratic 
equation  is  formed  by  substituting  (1)  into  (2).  The  quadratic  is 
solved  using  the  quadratic  formula.  Should  both  points  prove  to  be  ima¬ 
ginary,  then  the  ray  fails  to  intersect  the  circle. 

The  chord  length  l““  is  calculated  using  the  distance  formula. 
These  chord  lengths  are  formed  by  all  rays  ij  vhich  intersect  a  given 
parameter  circle  around  the  point  mn.  The  total  length  of  ray  segments 
dmn,  around  the  parameter  mn,  is  the  sum  of  all  these  chord  lengths. 


d“n  = 


(3) 


Every  ray  i j ,  vhich  passes  through  a  parameter  circle  mn,  vill  have 
a  residence  time  t®“  vithin  that  circle  as  calculated  from  the  observed 
travel  time  T.^  and  a  veighting  factor  v*° . 


mn 

Wij‘ 


(4) 


This  veighting  factor  depends  on  the  background  velocity  and  can  be 
vieved  as  the  ratio  of  the  expected  time  spent  inside  the  circle  to  the 
total  expected  time  of  the  ray.  For  a  constant  c  velocity  background: 


The  total  time  spent  abont  a  parameter  point  mn  is  the  sun  of  all 
the  residence  times.  The  average  slowness  at  the  point  mn  is  given  by: 


An  iterative  technique  is  applied  to  improve  upon  this  first  gness 
model.  First,  a  forward  model  is  needed  to  compote  an  estimated  travel 
time.  Using  straight  raypaths.  the  ray  is  divided  into  equal  segments 
Al.  The  slowness  along  the  entire  segment  length  is  assumed  to  equal 
the  slowness  at  the  midpoint  of  the  segment,  which  is  estimated  by  fit¬ 
ting  a  two-dimensional  parabola  to  the  surrounding  four  points  of  the 
slowness  matrii.  The  estimated  travel  time  T^  is  the  sum  of  the  slow¬ 


ness  along  the  raypath. 


(7) 


T  *  )  s(x,x)  Al. 

J  (xf.O) 


The  total  residual  travel  time  E, .  along  ray  ij  ia  given  by: 


E 


ij 


(8) 


This  residnal  error  time  ia  equally  distributed  along  the  length  of 
the  ray.  The  slowness  correction  Matrix  is  tbe  ratio  of  the  total  error 
tines  about  a  parameter  to  the  total  length  of  the  ray  segments. 


This  slowness  correction  matrix  is  added  to  the  previous  slowness 
matrix  to  obtain  an  updated  model.  Iteration  is  performed  by  computing 
a  new  correction  matrix,  based  on  the  current  slowness  model,  until  the 
mean  square  of  the  residual  times  reaches  a  steady  state. 


<  e. 


(10) 


2.i  Llaear  £(*)  Medina 


For  a  linear  c(z)  aedina,  the  raypaths  are  no  longer  straight,  bnt 
instead  are  curved.  The  local  radios  of  curvature  is  derived  using 
Snell 's  Lav. 


■M'*M  *  p  *  constant.  (11) 

c(z) 

Taking  the  derivative  of  the  ray  paraaeter,  and  noting  that  the 
derivative  of  any  constant  is  zero,  ve  find  that 


*  *  If  di  *  If 


(12) 


or 


0 


di 

c(z) 


-tin  i  ishl 

c*(z)  dz 


dz. 


Fron  the  geoaetry  shown  in  Fig.  3: 


dS  -  R di. 


and 


cos  i 


(13) 


(14) 


By  coahining  the  two  equations  (14) ,  ve  conclude  that 


di 


_ il  -  . 

R  cos  i 


(15) 


,\V  V  ,«  V  V  \*  *_■ 


Substituting  (15)  into  (13): 


Kc(z 


_ ii£_i  dc(l)  d 

)  c  *  ( z )  dz  ' 


1  .  tin i  jsdll  m  .  dc(z) 


R  c(z)  dz 


dz  ‘ 


So  for  «  linear  velocity  c(z)  *  cfl  +  «Dz,  the  radios  of  curvature 
becomes  a  constant,  which  implies  that  the  raypath  is  simply  the  arc 
along  some  circle;  that  is. 


I  *  ho¬ 


using  the  Wiechert-Herglotz  integral  (Aki  and  Richards.  1980),  the 
ray  parameter  can  be  determined  (See  appendix  A).  The  resnlt  is 


_ 2"o*i _ 

Vim  **.*  +  c  *  +  c*(z,)]  *  -  4c  *c*(z.) 
oio  j  o  j 


From  the  ray  parameter,  the  radius  of  curvature  becomes  known  using 
equation  (18).  Vith  two  points  on  this  circular  path  known,  namely  the 
shot  and  receiver  locations,  the  ray  follows  a  path  on  the  circle: 


(20) 


(x-xj*  ♦  <r  -  x  )*  «  R*. 
0  o 


Using  the  geoaetry  in  Fig.  4,  the  center  of  the  circle  can  be 
deterained: 


x  |  *  R  cos  (i+a) ; 


(21) 


i  J  *  I  sin  i. 


(22) 


Both  angles  are  calculated  froa  the  geoaetry.  The  resnlts  are 


.uf 


v  -  V 

2  R 


(23) 


tan (i + J)  «  “ 


(24) 


j 


Thus  for  a  linear  c(z)  nedina.  the  rays  are  arcs  of  a  circle.  The 
overall  length  of  the  ray  is: 


L 


ij 


-  Ra. 


(25) 


Over  the  tvo-diaensional  aedina,  a  grid  of  overlapping  circles  is 


superimposed.  A ny  ray  vill  inter  sect  e  circle  provided  that: 

| d  —  R |  <  r,  where  d*  ■  ( *B”*0 +  ^*n”*o^**  (26) 

The  arc  length  l"“  formed  by  the  ray  ij,  within  the  circle  an.  is 
computed  using  the  law  of  cosines  (Fig.  5). 

•••!  -  (27> 

1“  =  Rp.  (28) 

The  total  ray  distance  ^n.  about  a  given  parameter,  is  the  sub  of 
all  the  arc  lengths. 


d"n  - 


(29) 


Every  ray  ij,  which  passes  through  a  parameter  circle  an.  will  have 
a  residence  tiae  t"°  within  that  circle  as  calculated  from  the  observed 
travel  tiae  and  a  weighting  factor  w"°» 


wij- 


(30) 


The  weighting  factor  is  the  ratio  of  the  expected  tiae  spent  inside 


the  parameter  circle  to  the  overall  travel  tiae  for  the  entire  ray. 


boloid  to  the  surrounding  four  points  of  the  slowness  matrix.  An 
estimated  travel  tine  is  computed  by  sunning  the  slownesses  at  the  aid- 
points  of  the  segments. 


\  s(x,r)  Al. 

■*  A\ 


(xj.O) 


The  total  residual  time  along  the  ray  ij  is  given  by: 


-  V 


This  residual  error  time  is  distributed  equally  along  the  raypath. 
The  slowness  correction  aatriz  is  the  ratio  of  the  sua  of  the  local 
residual  tines  to  the  sun  of  the  arc  lengths  around  each  paraaeter: 


T  mn 

2  ij 
ii - , 

\lmn 
6  iJ 


jan 


This  slowness  correction  aatriz  is  added  to  the  slowness  aatriz  to 
obtain  a  revised  slowness  model.  Iteration  is  performed,  by  computing  a 
new  correction  matrix  based  on  the  updated  slowness  model,  until  the 


mean  square  of  the  residual  times  reaches  a  steady  state: 


IS 


(37) 


£(V  -  j><V*  <  - 

old  new 


1.  TESTING  TOE  ALGORITHMS 

Both  the  constant  and  linear  velocity  background  ART  algorithms 
were  tested  using  synthetic  data  sets.  For  a  given  slowness  nodel.  the 
travel  tines  were  generated  uaing  the  forward  modeling  scheme  in  either 
the  constant  or  linear  velocity  caae.  The  alowneas  was  integrated  along 
either  straight  or  circular  patha  as  dictated  by  the  general  velocity 
trend  of  the  input  velocity  model.  Vhere  the  velocity  increased 
linearly  with  z  inatead  of  z,  the  reciprocity  principle  of  interchanging 
the  shot  and  receiver  locations  was  applied.  Such  data  was  calculated 
by  interchanging  the  z  and  z  variables  and  using  the  technique  of  sec¬ 
tion  2.4  to  compute  the  synthetic  data. 

Clearly,  the  synthetic  data  does  not  correspond  ezactly  to  the  per¬ 
turbations  which  would  be  measured  in  the  real  case,  but  instead,  the 
travel  times  represent  a  first  order  approzimation  for  relatively  small 
amplitude  velocity  anomalies.  Therefore,  my  results  test  the  accuracy 
of  the  reconstruction  method  and  not  the  validity  of  linearizing  the 
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pro blea. 


4.  RESULTS 

ART  prodaces  the  best  tonogrsphical  resalts  when  the  velocity  func¬ 
tion  is  smoothly  varying.  Using  synthetic  data  sets,  accurately 
inverted  velocity  values  occurred  vithin  the  model,  but  the  velocity 
values  varied  from  the  synthetic  input  near  the  edges  of  the  model  (Fig. 
6,  7,  8,  9.  10.  11). 

The  greatest  deviation  occurred  near  the  origin  (Fig  6b.  7b,  8b. 

lib,  11c).  The  first  velocity  down  the  veil  is  erroneously  high  while 
the  first  velocity  along  the  surface  is  erroneously  low.  The  apparent 
cause  of  this  phenomenon  rests  within  the  raypaths  themselves  and  the 
way  in  which  the  ART  algorithm  handles  these  raypaths. 

The  first  parameter  down  the  well  reflects  the  slownesses  along  the 
rays  from  all  ten  shot  locations.  However,  the  shortest  raypath  length 
most  accurately  reflects  the  true  velocity  at  this  point,  while  the 
longer  rays  reflect  less  of  the  velocity  in  question  and  more  of  the 
velocities  further  away.  Because  the  ART  algorithm  equally  weights  all 
rays,  regardless  of  their  relative  lengths,  the  slowness  down  the  hole 
reflects  the  velocity  change  across  the  surface  of  the  model.  Appropri¬ 
ate  constraints  down  the  well  will  eliminate  this  problem,  since  the 
slowness  in  a  borehole  are  usually  known  from  the  sonic  log  (Fig  lid). 


Along  the  edges,  the  velocity  velnes  fluctuated  from  the  synthetic 
■odel.  These  deviations  resulted  from  two  factors.  First,  the  errone¬ 
ous  origin  velocities  greatly  influence  the  surface  and  borehole  edges. 
Second,  the  diagonal  edge  parameters  are  poorly  constrained. 

Unconstrained  paraaeters  are  defined  by  only  one  ray.  All  parame¬ 
ters  within  my  model  had  at  least  three  rays  defining  each  parameter. 
However,  along  the  diagonal  edge,  these  three  rays  are  nearly  identical 
and  as  such,  they  are  defined  by  nearly  the  same  parameters.  Therefore, 
these  rays  fail  to  be  completely  independent  of  each  other.  Thus,  the 
parameters  along  the  diagonal  edge  are  poorly  constrained.  In  uncon¬ 
strained  cases,  ART  tends  to  average  the  velocities  equally  along  the 
raypath  resulting  in  resolution  problems. 

By  using  a  delta  function  as  the  input  velocity  structure  (Fig  9), 
we  can  better  illustrate  this  resolution  problem.  Resolution  along  the 
raypaths  is  poor,  but  between  the  rays,  the  resolution  increases.  This 
accounts  for  the  skewness  of  the  inverted  delta  function  into  an  ellipt¬ 
ical  form  elongated  along  the  dominant  raypath  direction.  Still,  the 
location  of  the  slowness  delta  point  is  accurately  defined,  though  not 
its  true  velocity  value. 

This  poor  resolution  along  the  raypaths  contributes  to  uniqueness 
problems  (Fig  10).  In  general,  an  infinite  number  of  models  can  produce 
the  same  travel  time  data.  ART  converges  to  the  model  which  is  the 
closest  to  the  assumed  background  velocity  and  not  necessarily  a 


representative  of  the  trne  medium. 

Uniqueness  requires  that  all  the  slowness  parameters  are  over  con¬ 
strained;  that  is  more  than  one  ray  intersects  each  parameter  circle. 
Furthermore,  these  rays  must  be  independent.  A  close  examination  of  my 
model  reveals  that  at  least  three  rays  intersect  each  parameter;  how¬ 
ever,  these  rays  travel  approximately  along  the  came  direction  and  thus, 
these  rays  reflect  many  of  the  same  slowness  parameters.  Clearly,  these 
rays  are  not  independent  of  each  other. 

The  back-projected  algorithm  along  curved  rays  was  tested  nsing 
synthetic  data  generated  on  linear  trending,  velocity  models.  In  each 
case,  the  curved  ART  algorithm  accurately  reconstructed  the  synthetic 
models  (Fig.  12,  13,  14). 

5.  PIERRE  SHALE  EXPERIMENT 

The  multiple  offset  VSP,  which  was  taken  by  the  CSM  Exploration 
Research  Laboratory  of  the  Geophysics  Department,  provided  a  real  data 
test  of  the  algorithms.  Field  work  was  done,  in  November  1980,  at  the 
CSM  test  site  in  northeastern  Colorado,  about  five  miles  south  of  Brush 
Colorado  (Southwest  1/4,  Section  28,  Township  3N,  Range  55V).  Sixta 
(1982)  describes  the  data  aquisition  and  data  processing  of  this  VSP 
experiment. 

This  site  was  chosen  for  its  simple  geologic  setting.  A  40  foot 


thick  layer  of  eolian  sand  rests  on  the  surface.  Beneath  this  sand 
layer,  a  40  foot  thick  layer  of  clay  grades  into  the  Pierre  Shale  below. 
These  thicknesses  vary  laterally  away  from  the  borehole.  Because  the 
npper  1000  feet  of  the  Pierre  Shale  is  well  noted  for  its  homogeneous 
nature,  refracted  events  should  not  appear  on  the  records  and  therefore, 
the  direct  ray  will  be  the  first  arrival. 

The  field  layout  is  sketched  in  Fig.  IS.  Geophones  were  spaced 
down  the  cemented  well  every  100  feet  within  the  interval  between  400 
and  1000  feet.  Dynamite,  placed  in  shot  holes  at  a  depth  of  100  feet 
deep,  were  shot  at  200  foot  intervals  across  the  surface  to  a  distance 
of  1400  feet  from  the  borehole.  The  travel  times  were  picked  at  the 
first  trough  on  the  records  (Fig.  16)  and  are  tabulated  in  the  table 
below. 


Table  1:  Travel  times  in  msec. 

Horizontal  offset  in  ft. 


Depth 
in  ft. 

200 

400 

600 

800 

1000 

1200 

1400 

400 

50.25 

62.00 

82.75 

106.25 

133.50 

160.00 

187.00 

500 

64.50 

74.00 

91.75 

113.00 

138.50 

163.50 

189.75 

600 

78.75 

86.50 

101.75 

121.00 

144.50 

168.00 

193.50 

700 

92.25 

99.00 

112.25 

129.75 

151.50 

174.00 

198.25 

800 

106.00 

111.75 

123 .50 

139.25 

159.50 

180.75 

203.50 

900 

119.25 

124.50 

135.00 

149.00 

168.00 

188.00 

210.00 

1000 

133 .50 

138.00 

147.00 

160.00 

177.75 

196.50 

217.00 

Picking  the  arrival  times  at  the  first  trough  results  in  a  lag  time 
of  about  3  msec,  from  the  true  onset  of  the  waveform.  This  error  in  the 
travel  time  corresponds  to  a  similar  error  in  the  velocity  value.  The 


travel  tine  is  the  integral  of  the  slowness  along  the  raypath. 


t  =  f  —  dL,  or  t  -  k  (38) 

J  c  c 

Differentiating  (38),  we  obtain 

dt  *  -L^f .  (39) 

c* 

Dividing  (38)  into  (39),  yields  the  result 

41  «  -L^t.  (40) 

t  c 

By  setting  dt=3  nsec,  and  noting  that  t  ranges  from  50.25  to  217.0 
■sec.,  the  true  velocity  values  range  fron  6.0  to  1.3%  faster  than  those 
which  are  shown  in  the  reconstructed  fields. 

The  Pierre  Shale  is  well  noted  for  its  honogeneous  nature;  however, 
our  reconstructions  fail  to  support  this  characteristic.  Imaging  was 
performed  both  without  paraaeter  constraints  and  with  constraining  the 
slowness  parameters  down  the  borehole  to  within  five  percent  of  the 
measured  slowness  values  off  the  sonic  log. 

Both  the  constrained  (Fig.  18a,  18b)  and  the  unconstrained  (Fig. 
17a,  17b)  reconstructions  reveal  a  questionable  high  velocity  cell  near 
the  surface.  Although  the  sonic  log  for  the  well  records  a  velocity  of 
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6250  ft/sec  at  a  depth  of  100  feet,  the  reconstructed  high,  only  600 
feet  away,  reaches  11,900  ft/sec  for  the  unconstrained  case  and  12,000 
ft/sec  for  the  constrained  case.  However,  these  velocity  highs  con¬ 
sisted  of  a  single,  poorly-constrained  parameter  along  the  surface, 
which  clearly  introduces  uncertainty  into  this  velocity  value. 

By  eliminating  any  poorly  constrained  parameters  along  the  surface, 
the  imaging  process  is  improved  (Fig.  19,  20).  These  errant  velocity 
values  were  ignored  only  when  the  reconstructed  velocity  field  was  con¬ 
toured,  and  thus,  they  were  used  throughout  the  reconstruction  algo¬ 
rithm.  A  similar  step  was  applied  to  the  reconstructions  using  circular 
rays  based  on  the  linear  velocity  trend  exhibited  by  the  sonic  log  (Fig. 
21,  22).  Each  reconstruction  shows  a  high  velocity  ridge  of  8,000 
ft/sec.  hovering  between  a  horizontal  offset  of  200  to  900  feet  at  a 
depth  around  200  to  300  feet.  Additionally,  a  low  velocity  cell  exists 
at  100  foot  offset  and  a  depth  of  500  to  800  feet.  This  low  cell  fluc¬ 
tuates  between  velocities  of  7,000  to  7,200  ft/sec.  depending  upon  the 
reconstruction  conditions. 

To  test  the  validity  of  these  reconstructions,  the  observed  travel 
times  were  compared  to  theoretical  data  from  a  linearly  increasing 
medium.  This  synthetic  model,  based  on  sonic  log  information,  initially 
started  at  a  velocity  of  7,000  ft/sec.  and  increased  to  8,000  ft/sec.  at 
the  total  depth  of  the  well.  For  comparision,  both  the  synthetic  and 
the  Pierre  Shale  VSP  travel  times  were  plotted  for  each  geophone  poai- 
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tion  (Fig.  23).  Where  the  synthetic  carve  arrives  after  the  observed 
tines,  a  velocity  increase  is  required.  Where  the  observed  tiaes  lag 
behind  the  synthetic,  lover  velocities  are  expected.  Relative  to  the 
synthetic  aodel.  ve  should  expect  higher  velocities  at  the  large  offset 
distances,  but  lover  velocities  at  snail  offsets  and  large  depth. 


Because  surface  velocities  are  typically  very  slov,  the  high  sur¬ 
face  velocity  is  contrary  to  expectations.  Three  feasible  explanations 
arise  for  this  velocity  anoaaly.  Either  the  velocity  structure  exists, 
or  the  data  or  the  collection  of  the  data  is  questionable,  or  the  basic 
theory  is  inadequate  for  the  earth  conditions  exhibited  by  the  Pierre 
Shale.  Poorly  placed  shot  locations  could  result  in  an  artificial  high 
at  the  surface,  if  the  errant  shot  locations  resulted  in  a  shorter  ray- 
path. 


A  aore  feasible  explanation  arises  from  the  strong  anisotropy  exhi¬ 
bited  by  the  Pierre  Shale.  After  using  the  sane  data  collected  near 
Brush,  Colorado,  White,  Martineau-Nicoletis,  and  Monash  (1983)  concluded 
that  the  horizontal  velocity  component  vss  10  to  20%  faster  than  the 
vertical  velocity  component.  Qualitatively,  the  rays  froa  the  large 
offset  sources  possess  mostly  a  horizontal  coaponent  and  thus,  should 
exhibit  a  greater  than  expected  velocity.  Similarly,  the  deeper  travel¬ 
ing  rays  are  mostly  vertical  and  vould  travel  at  the  slower  velocity. 
Therefore,  the  anisotropy  of  the  Pierre  Shale  explains  the  high  surface 
velocity  and  the  low  velocity  cell  at  depth  of  the  reconstructed  velo- 


city  fields. 


£.  CONCLUSIONS 

Algebraic  reconstruction  technique  (ART)  concentrates  on  the  pro¬ 
duction  of  a  reconstructed  field  whose  projected  data  (travel  tines) 
agree  with  the  observed  data.  This  reconstructed  field  is  modified  by 
altering  the  data  for  each  ray  such  that  when  this  data  is  back- 
projected,  the  new  inage  agrees  with  the  original  data. 

ART  possesses  aany  inherent  advantages.  The  flexibility  of  the 
algorithm  easily  allows  for  the  incorporation  of  constraints.  For  exam¬ 
ple,  sonic  log  information  predetermines  the  velocity  parameters  down 
the  borehole.  This  flexibility  also  permits  the  application  of  ART  to 
any  field  setup  of  sources  and  receivers. 

With  noiseless  data,  ART  picks  a  projection  which  agrees  with  the 
observed  data.  It  is  not  in  ART’s  nature  to  introduce  spurious  images 
on  good  data  sets.  However,  the  resolution  of  the  image  is  affected  by 
the  orientation  of  the  anomaly  relative  to  the  direction  of  the  rays. 
While  good  resolving  power  occurs  perpendicular  to  the  rays,  ART  tends 
to  smooth  the  image  along  the  direction  of  the  rays. 

With  noisy  data,  ART  can  produce  single  point  anomalies  of  highly 
unrealistic  velocities  among  the  poorly-constrained  parameters  along  the 
edges  of  the  model.  ART  will  produce  a  reasonable  reconstructed  velo- 


city  field  when  constraints  are  incorporated  from  the  well  log  informa¬ 
tion  or  by  simply  ignoring  any  poorly-constrained  parameters  within  the 
image.  Lastly,  data,  which  was  collected  over  a  strongly  anisotropic 
medium,  will  result  in  poor  velocity  reconstructions,  because  such  data 
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APPENDII  A:  RAT  PARAMETER  DERIVATION 


The  ray  parameter  within  a  linear  c(z)  aediua  is  determined  froa 


the  Wiechert-Herglotz  integral  (Aki  and  Richards.  1980). 


'i  *  P  Jl  ^ 

1  JQ  Vn*  -  p* 


(A-l) 


where  ij  *  **  •!<>*»***» 


is  the  source  offset. 


Zj  is  the  receiver  depth. 


The  aedinm  is  assnaed  to  possess  a  linear  velocity  structure. 


c(z)  =  c  +  ■  z;  c..  “  c(0),  the  surface  velocity.  (A-2) 

o  o  o 


Substituting  (A-2)  into  (A-l),  we  obtain 


/J  <ct,+ll0*)d* 

i.  *  p  I  . 

1  0  Vl  -  p»  ( C  +  B_zj* 


(A-3) 


o  o 


Hie  integration  aay  be  carried  out  explicitly.  The  result  is 


z“2j 

x.  *  — Vl  -  p*c*(z)  I  , 

*  ■  ©  1 _ A 


(A-4) 
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«  *  .  '  .  w  .  •  a"  .  •  •  *  •  a' 


■i  ■  i1;  [v*  -  p***<v  -  vi  -  ■>*%*] 


(A-5) 


Squaring  both  sides,  we  find 


(xi*,“o)*  "  (1-P*c0»)  +  (l-p>c»(ij))  -  2  V(l-p*co»)  (l-p»c»(Xj))  .  (A-6) 


or 


(xiP“o)*  "  2  +  P*(co*  +  c*(lj)>  “  -2  V(1  -  P*c0*)  (1  -  p*c»(Zj)).  (A-7) 

Again  squaring  both  sides,  we  find 

P4  I®  *x.  *  +  e  *  +  e*(z  ,)]  ~  A  «  *p*l«  *x.  *  +  c  *  +  c*(  i.)]  +  4  *  (A-8) 

oio  j  o'oio  j 

4(l-p»co*)  (l-p»c»(zj)). 

or 

P4[  [■0*xi*  +  c0*  +  c*(Zj)]*  -  4ce*c*(Zj)]  -  4p*m0*zi*  -0.  (A-9) 

Solving  for  p,  we  obtain 

_ 4ao>xi> _ 

[®_*x .  *  +  cft*  +  c*(  z4)  ]  *  -  4c  ‘cMz.)’ 


p*  *  0. 


(A-10) 


However,  the  rey  parameter  must  be  positive  for  a  linearly  increas¬ 


ing  Medina.  Thu*,  we  conclnde  that 


APPENDIX  B:  TRAVEL  TIME  CALCULATION 


The  travel  time  is  a  linear  c(i)  medium  is  calculated  by  the 


Wiechert-Herglotz  inversion  method  (Aki  and  Richards,  1980). 


f  ' - - . 

0  c(z)  Vl  -  p*c*(z) 


We  introduce  the  new  variable  of  integration  0: 


pc(z)  =  sinO, 


^1  -  p*c*(z)  *  cos  9. 


Differentiating  (B-3),  we  obtain 


pmQ  dz  *  cos  9  d9. 


Substituting  (B-2)  and  (B-4)  into  (B-l)  yields: 


T  J_  f  d9 

ij  =  m  J  s in  9 ' 


(B-l) 


(R-2) 


(B-3) 


(B-4) 


(B-5) 
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Travel  tinea  calculated  using  curved  rays  based  After  7  iterations  using  straight  rays 
on  a  background  velocity  -  v(x)  **  5000  ♦  .01(*.). 


background  velocity  -  v(x)  “  5000  +  .l(x^). 

Figure  7.  Velocity  increasing  with  x 


Figure  9.  Model  of  ■  low  Telocity  anouly 


Pignr«  11.  Tvo-lifertd  aodel 


Reconstructed  Velocity  Field 


Figure  11.  Two-layered  model 


Figure  12.  Dipping  model. 


Figure  13.  High  velocity  anomaly  within  a  linear  trend. 
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Figure  15.  Field  experiment  layout. 


After  1  iteration  using  straight  rays.  Velocity  parameters  were 
constrained  down  the  well  (z-axis). 


Reconstructed  Velocity  Field 


After  5  iterations  using  straight  rays.  Poorly-constrained  paraaeters 
along  the  surface  were  eliminated  from  the  image. 


Reconstructed  Velocity  Field 


Figure  20b.  Acoustic  tonography  from  the  Pierre  Shale  VSP  data 


Reconstructed  Velocity  Field 


Figure  22a.  Acoustic  tomography  from  the  Pierre  Shale  VSP  data. 
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Figure  23.  Travel  time  curves. 
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An  algebraic  reconstruction  technique  (ART)  is  described  for  the 
seismic  tomogrsphy  of  velocities  from  the  travel  times  of  a  multiple 
offset  vertical  seismic  profile.  ART  concentrates  on  the  production  of 
a  reconstructed  field  whose  projected  data  (travel  times)  agree  with  the 
observed  date.  This  reconstructed  field  is  modified  by  altering  the  data 
for  each  ray  such  that  when  this  data  is  back-projected,  the  new  image 
agrees  with  the  original  data.  Because  the  paths  of  the  rays  must  be 
known  to  calculate  the  expected  travel  times,  the  problem  is  linearised 
by  using  raypath  approximations  as  determined  from  either  a  constant  or 
a  linear  c(z)  velocity  medium. 


Imaging  of  synthetic  data  revealed  that  the  orientation  of  the  ano¬ 
maly  affects  both  the  rate  of  convergence  and  the  resolution  of  the 
reconstructed  field.  Some  smoothing  of  the  velocity  anomalies  occurred 
along  the  direction  of  the  rays. 


Noisy  data  sets  developed  problems  in  the  reconstructed  velocity 
field.  Huge  single  point  anomalies  appeared  along  the  model's  edge  in 
the  reconstructed  image.  Elimination  of  these  errant  anomalies  was 
necessary  to  obtain  a  reasonable  velocity  reconstruction.  Because  iso¬ 
tropy  was  assumed,  the  algorithm  poorly  reconstructed  data  which  was 
collected  over  the  strongly  anisotropic  Pierre  Shale. 
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